Statistics Fundamentals: interactive
PLS 206 — Applied Multivariate Modeling
Apps for Week 1, in the order the lecture covers them. They run in the browser, so no R installation is needed. Tick “Show R code” in any app to get code that reproduces the current view in RStudio.
Distribution Explorer
Choose a distribution family, adjust its parameters, apply a transform, and see the histogram, boxplot, and Q–Q plot update together. The normal overlay shows whether the data (or its transform) is approximately normal.
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 920
library(shiny)
ui <- fluidPage(
sidebarLayout(
sidebarPanel(
selectInput("dist", "Distribution",
choices = c("Normal" = "norm",
"Uniform" = "unif",
"Exponential" = "exp",
"Log-normal" = "lnorm",
"t (heavy tails)" = "t")),
conditionalPanel("input.dist == 'norm'",
sliderInput("norm_mean", "Mean (μ)", -10, 10, 0, 0.5),
sliderInput("norm_sd", "Std Dev (σ)", 0.1, 5, 1, 0.1)
),
conditionalPanel("input.dist == 'unif'",
sliderInput("unif_min", "Min", -10, 10, 0, 0.5),
sliderInput("unif_max", "Max", 0, 20, 5, 0.5)
),
conditionalPanel("input.dist == 'exp'",
sliderInput("exp_rate", "Rate (λ = 1/mean)", 0.1, 5, 1, 0.1)
),
conditionalPanel("input.dist == 'lnorm'",
sliderInput("ln_meanlog", "Mean of log(x)", -2, 3, 1, 0.1),
sliderInput("ln_sdlog", "SD of log(x)", 0.1, 2, 0.5, 0.1)
),
conditionalPanel("input.dist == 't'",
sliderInput("t_df", "Degrees of freedom", 1, 30, 3, 1)
),
sliderInput("n", "Sample size (n)", 20, 1000, 200, 20),
hr(),
selectInput("transform", "Transform x",
choices = c("x (no transform)" = "none",
"log(x)" = "log",
"sqrt(x)" = "sqrt",
"x^2" = "sq",
"1/x" = "inv")),
hr(),
actionButton("resample", "New sample"),
checkboxInput("show_code", "Show R code", FALSE),
helpText("Try: Log-normal + log(x). The Q-Q plot becomes a straight line."),
helpText("Then Exponential + log(x): the skew flips to the left instead of going away.")
),
mainPanel(
plotOutput("distPlot", height = "530px"),
verbatimTextOutput("stats"),
conditionalPanel("input.show_code == true",
verbatimTextOutput("code_out")
)
)
)
)
server <- function(input, output, session) {
seed <- reactive(206 + input$resample)
raw_data <- reactive({
n <- input$n
set.seed(seed())
switch(input$dist,
norm = rnorm(n, input$norm_mean, input$norm_sd),
unif = runif(n, min(input$unif_min, input$unif_max - 0.1),
max(input$unif_max, input$unif_min + 0.1)),
exp = rexp(n, input$exp_rate),
lnorm = rlnorm(n, input$ln_meanlog, input$ln_sdlog),
t = rt(n, input$t_df)
)
})
plot_data <- reactive({
x <- raw_data()
switch(input$transform,
none = x,
log = { v <- x[x > 0]; log(v) },
sqrt = { v <- x[x >= 0]; sqrt(v) },
sq = x^2,
inv = { v <- x[x != 0]; 1/v }
)
})
x_label <- reactive({
switch(input$transform,
none = "x", log = "log(x)", sqrt = "sqrt(x)", sq = "x\u00b2", inv = "1/x"
)
})
output$distPlot <- renderPlot({
x <- plot_data()
lbl <- x_label()
if (length(x) < 5) {
plot.new()
text(0.5, 0.5,
"Too few valid values after transform.\nTry a different distribution or transform.",
cex = 1.2)
return(invisible(NULL))
}
layout(matrix(1:3, ncol = 1), heights = c(1, 2.5, 2.5))
# 1. Boxplot
par(mar = c(1, 4, 2, 1))
boxplot(x, horizontal = TRUE,
col = adjustcolor("steelblue", 0.45), border = "#2c5f8a",
xaxt = "n",
main = if (input$transform != "none")
paste0(input$dist, " \u2192 transform: ", lbl)
else input$dist)
# 2. Histogram + normal overlay
par(mar = c(3, 4, 0.5, 1))
h <- hist(x, plot = FALSE, breaks = 30)
ymax <- max(h$density) * 1.35
hist(x, freq = FALSE, breaks = 30,
col = adjustcolor("red", 0.45), border = "white",
ylim = c(0, ymax), xlab = lbl, ylab = "Density", main = "")
mu <- mean(x); sg <- sd(x)
xs <- seq(min(x) - 0.1 * diff(range(x)),
max(x) + 0.1 * diff(range(x)), length.out = 300)
lines(xs, dnorm(xs, mu, sg), col = "#002855", lwd = 2)
abline(v = mu, col = "orange", lty = 2, lwd = 1.8)
legend("topright", bty = "n", cex = 0.8,
legend = c(sprintf("Normal(\u03bc\u0302=%.2f, \u03c3\u0302=%.2f)", mu, sg),
"Mean"),
col = c("#002855", "orange"), lty = c(1, 2), lwd = 2)
# 3. Q-Q plot
par(mar = c(4, 4, 0.5, 1))
qqnorm(x, pch = 16, cex = 0.5,
col = adjustcolor("steelblue4", 0.6),
main = "", xlab = "Theoretical normal quantiles",
ylab = paste0("Sample quantiles (", lbl, ")"))
qqline(x, col = "orange", lwd = 2)
})
output$stats <- renderPrint({
x <- plot_data()
raw <- raw_data()
lbl <- x_label()
if (length(x) < 5) {
cat("Too few valid values after transform.\n"); return(invisible(NULL))
}
n_raw <- length(raw); n_v <- length(x)
if (n_v < n_raw)
cat(sprintf("Note: %d/%d values filtered by transform (domain restriction).\n\n",
n_raw - n_v, n_raw))
q <- quantile(x, c(0.25, 0.5, 0.75))
cat(sprintf("n = %d\n", n_v))
cat(sprintf("Mean = %8.3f | Median = %.3f\n", mean(x), q[2]))
cat(sprintf("SD = %8.3f | IQR = %.3f\n", sd(x), q[3] - q[1]))
cat(sprintf("Q1 = %8.3f | Q3 = %.3f\n", q[1], q[3]))
cat(sprintf("Min = %8.3f | Max = %.3f\n\n", min(x), max(x)))
sw <- tryCatch(shapiro.test(if (n_v > 5000) sample(x, 5000) else x),
error = function(e) NULL)
if (!is.null(sw)) {
cat(sprintf("Shapiro-Wilk: W = %.4f, p = %.4f\n", sw$statistic, sw$p.value))
cat(if (sw$p.value < 0.05)
" \u2192 Evidence against normality (p < 0.05)\n"
else
" \u2192 Consistent with normality (p \u2265 0.05)\n")
}
})
output$code_out <- renderPrint({
lbl <- x_label()
n <- input$n
gen <- switch(input$dist,
norm = sprintf("x <- rnorm(%d, mean = %.1f, sd = %.1f)",
n, input$norm_mean, input$norm_sd),
unif = sprintf("x <- runif(%d, min = %.1f, max = %.1f)",
n, min(input$unif_min, input$unif_max - 0.1),
max(input$unif_max, input$unif_min + 0.1)),
exp = sprintf("x <- rexp(%d, rate = %.1f)", n, input$exp_rate),
lnorm = sprintf("x <- rlnorm(%d, meanlog = %.1f, sdlog = %.1f)",
n, input$ln_meanlog, input$ln_sdlog),
t = sprintf("x <- rt(%d, df = %d)", n, input$t_df)
)
trans <- switch(input$transform,
none = "", log = "x <- log(x[x > 0])", sqrt = "x <- sqrt(x[x >= 0])",
sq = "x <- x^2", inv = "x <- 1 / x[x != 0]"
)
cat("# Minimal R code to reproduce this view:\n\n")
cat(sprintf("set.seed(%d)\n", seed()))
cat(gen, "\n")
if (nchar(trans) > 0) cat(trans, "\n")
cat("\npar(mfrow = c(3, 1))\n")
cat(sprintf("boxplot(x, horizontal = TRUE, xlab = \"%s\")\n", lbl))
cat(sprintf("hist(x, freq = FALSE, breaks = 30, xlab = \"%s\")\n", lbl))
cat(sprintf("qqnorm(x); qqline(x)\n\n"))
cat("summary(x)\nshapiro.test(x)\n")
})
}
shinyApp(ui, server)SD vs SE: Many Experiments
One experiment gives us one sample of plants and one mean. Here we repeat the experiment many times: the SD (left) describes how plants vary, the SE (right) describes how the mean would vary across repeated experiments. Raise n and watch only the right panel narrow, as SE = SD / sqrt(n).
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 760
library(shiny)
pops <- list(
norm = list(
gen = function(N) rnorm(N, mean = 82, sd = 16),
gen_txt = "rnorm(%d, mean = 82, sd = 16)",
mu = 82,
sigma = 16,
mu_txt = "82",
sig_txt = "16",
lim_txt = "qnorm(c(0.001, 0.999), mean = 82, sd = 16)",
lim = function() qnorm(c(0.001, 0.999), mean = 82, sd = 16),
xlab = "Height (cm)"
),
lnorm = list(
gen = function(N) rlnorm(N, meanlog = log(20), sdlog = 0.7),
gen_txt = "rlnorm(%d, meanlog = log(20), sdlog = 0.7)",
mu = exp(log(20) + 0.49 / 2),
sigma = sqrt((exp(0.49) - 1) * exp(2 * log(20) + 0.49)),
mu_txt = "exp(log(20) + 0.7^2 / 2)",
sig_txt = "sqrt((exp(0.7^2) - 1) * exp(2 * log(20) + 0.7^2))",
lim_txt = "qlnorm(c(0.001, 0.999), meanlog = log(20), sdlog = 0.7)",
lim = function() qlnorm(c(0.001, 0.999), meanlog = log(20), sdlog = 0.7),
xlab = "Biomass (g)"
),
exp = list(
gen = function(N) rexp(N, rate = 1/30),
gen_txt = "rexp(%d, rate = 1/30)",
mu = 30,
sigma = 30,
mu_txt = "30",
sig_txt = "30",
lim_txt = "qexp(c(0.001, 0.999), rate = 1/30)",
lim = function() qexp(c(0.001, 0.999), rate = 1/30),
xlab = "Wait (days)"
)
)
ui <- fluidPage(
sidebarLayout(
sidebarPanel(
selectInput("pop", "Population",
choices = c("Normal: height ~ N(82, 16^2)" = "norm",
"Lognormal biomass (meanlog = log(20), sdlog = 0.7)" = "lnorm",
"Exponential waits (mean 30)" = "exp")),
sliderInput("n", "Plants per experiment (n)", min = 2, max = 100, value = 12, step = 1),
sliderInput("k", "Number of experiments", min = 100, max = 2000, value = 1000, step = 100),
actionButton("resample", "New sample"),
hr(),
checkboxInput("show_code", "Show R code", FALSE),
helpText("Raise n: the left panel stays as wide, the right panel narrows."),
helpText("Pick Exponential with n = 2, then n = 30: the means become bell-shaped (CLT).")
),
mainPanel(
plotOutput("sePlot", height = "400px"),
verbatimTextOutput("stats"),
conditionalPanel("input.show_code == true",
verbatimTextOutput("code_out")
)
)
)
)
server <- function(input, output, session) {
seed <- reactiveVal(206)
observeEvent(input$resample, seed(seed() + 1))
sim <- reactive({
p <- pops[[input$pop]]
n <- input$n; k <- input$k
set.seed(seed())
m <- matrix(p$gen(k * n), nrow = k)
x <- m[1, ]
means <- rowMeans(m)
list(p = p, n = n, k = k, x = x, means = means,
xlim = range(x, means, p$lim()))
})
output$sePlot <- renderPlot({
s <- sim(); p <- s$p; x <- s$x; means <- s$means
par(mfrow = c(1, 2))
# LEFT: one experiment
hist(x, col = "steelblue", border = "white", xlim = s$xlim,
main = sprintf("One experiment (n = %d)", s$n), xlab = p$xlab)
abline(v = mean(x), col = "#002855", lwd = 3)
abline(v = p$mu, col = "orange", lwd = 3, lty = 2)
y1 <- par("usr")[4] * 0.9
arrows(mean(x) - sd(x), y1, mean(x) + sd(x), y1,
code = 3, col = "#FF0000", lwd = 2, length = 0.08)
text(mean(x), y1, "± 1 SD", pos = 3, col = "#FF0000", xpd = TRUE)
# RIGHT: means of all experiments
h2 <- hist(means, breaks = 30, plot = FALSE)
ymax <- max(h2$density, dnorm(mean(means), mean(means), sd(means))) * 1.15
hist(means, breaks = 30, freq = FALSE, col = "steelblue", border = "white",
xlim = s$xlim, ylim = c(0, ymax),
main = sprintf("Means of all %d experiments", s$k), xlab = paste("Mean", tolower(p$xlab)))
abline(v = p$mu, col = "orange", lwd = 3, lty = 2)
curve(dnorm(x, mean(means), sd(means)), add = TRUE, col = "#002855", lwd = 2, n = 400)
y2 <- ymax * 0.93
arrows(mean(means) - sd(means), y2, mean(means) + sd(means), y2,
code = 3, col = "#00B050", lwd = 2, length = 0.05)
text(mean(means), y2, "± 1 SE", pos = 3, col = "#00B050", xpd = TRUE)
})
output$stats <- renderPrint({
s <- sim(); x <- s$x
cat(sprintf("True population SD (sigma): %8.3f\n", s$p$sigma))
cat(sprintf("SD of experiment 1, sd(x): %8.3f\n", sd(x)))
cat(sprintf("SE from that one sample: %8.3f (sd(x) / sqrt(n))\n", sd(x) / sqrt(s$n)))
cat(sprintf("SD of the means (simulated SE): %8.3f (sd(means), %d experiments)\n", sd(s$means), s$k))
cat(sprintf("Theory, sigma / sqrt(n): %8.3f (n = %d)\n", s$p$sigma / sqrt(s$n), s$n))
cat("\nSE = SD / sqrt(n)\n")
})
output$code_out <- renderPrint({
s <- sim(); p <- s$p; n <- s$n; k <- s$k
cat("# R code to reproduce this view\n\n")
cat(sprintf("set.seed(%d)\n", seed()))
cat(sprintf("n <- %d # plants per experiment\n", n))
cat(sprintf("k <- %d # number of experiments\n", k))
cat("# one call makes all k experiments: each row is one experiment\n")
cat(sprintf("m <- matrix(%s, nrow = k)\n", sprintf(p$gen_txt, k * n)))
cat("x <- m[1, ] # experiment 1\n")
cat("means <- rowMeans(m) # the mean of every experiment\n\n")
cat(sprintf("mu <- %s # true mean\n", p$mu_txt))
cat(sprintf("sigma <- %s # true SD\n\n", p$sig_txt))
cat("sigma # true population SD\n")
cat("sd(x) # SD of experiment 1\n")
cat("sd(x) / sqrt(n) # SE estimated from that one sample\n")
cat("sd(means) # SD of the means = simulated SE\n")
cat("sigma / sqrt(n) # SE from theory\n\n")
cat(sprintf("xlim <- range(x, means, %s)\n", p$lim_txt))
cat("par(mfrow = c(1, 2))\n")
cat(sprintf("hist(x, col = \"steelblue\", border = \"white\", xlim = xlim,\n main = \"One experiment (n = %d)\", xlab = \"%s\")\n", n, p$xlab))
cat("abline(v = mean(x), col = \"#002855\", lwd = 3)\n")
cat("abline(v = mu, col = \"orange\", lwd = 3, lty = 2)\n")
cat("y1 <- par(\"usr\")[4] * 0.9\n")
cat("arrows(mean(x) - sd(x), y1, mean(x) + sd(x), y1, code = 3, col = \"#FF0000\", lwd = 2)\n\n")
cat("h2 <- hist(means, breaks = 30, plot = FALSE)\n")
cat("ymax <- max(h2$density, dnorm(mean(means), mean(means), sd(means))) * 1.15\n")
cat(sprintf("hist(means, breaks = 30, freq = FALSE, col = \"steelblue\", border = \"white\",\n xlim = xlim, ylim = c(0, ymax), main = \"Means of all %d experiments\")\n", k))
cat("abline(v = mu, col = \"orange\", lwd = 3, lty = 2)\n")
cat("curve(dnorm(x, mean(means), sd(means)), add = TRUE, col = \"#002855\", lwd = 2, n = 400)\n")
cat("y2 <- ymax * 0.93\n")
cat("arrows(mean(means) - sd(means), y2, mean(means) + sd(means), y2, code = 3, col = \"#00B050\", lwd = 2)\n")
})
}
shinyApp(ui, server)100 Experiments, 100 Intervals
Each row is one experiment: we measure n wheat plants from a field where height ~ N(82, 16^2) and build a confidence interval for the mean. Click “New sample” a few times and watch the running total settle near the confidence level: the 95% describes the method, while any one interval either catches 82 cm or misses it (Week 1 script, section 12).
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 900
library(shiny)
ui <- fluidPage(
sidebarLayout(
sidebarPanel(
sliderInput("n", "Plants per experiment (n)", min = 3, max = 50, value = 12, step = 1),
selectInput("level", "Confidence level",
choices = c("80%" = "0.80", "90%" = "0.90",
"95%" = "0.95", "99%" = "0.99"),
selected = "0.95"),
radioButtons("mult", "Multiplier",
choices = c("t: qt(1 - a/2, n - 1)" = "t",
"z: qnorm(1 - a/2) (wrong for small n)" = "z"),
selected = "t"),
helpText("The population is fixed: plant height ~ N(82, 16^2), so the true mean is 82 cm."),
actionButton("resample", "New sample"),
hr(),
checkboxInput("show_code", "Show R code", FALSE),
helpText("Switch to z with n = 3: far fewer than 95% cover."),
helpText("Raise the confidence level: intervals get wider and fewer miss.")
),
mainPanel(
plotOutput("ciPlot", height = "560px"),
verbatimTextOutput("stats"),
conditionalPanel("input.show_code == true",
verbatimTextOutput("code_out")
)
)
)
)
server <- function(input, output, session) {
mu <- 82
sigma <- 16
seed <- reactive(206 + input$resample)
sim <- reactive({
n <- input$n
level <- as.numeric(input$level)
a <- 1 - level
set.seed(seed())
m <- matrix(rnorm(100 * n, mu, sigma), nrow = 100)
xbar <- rowMeans(m)
se <- apply(m, 1, sd) / sqrt(n)
mult <- if (input$mult == "t") qt(1 - a / 2, n - 1) else qnorm(1 - a / 2)
lower <- xbar - mult * se
upper <- xbar + mult * se
covers <- lower <= mu & upper >= mu
list(n = n, level = level, a = a, m = m, xbar = xbar, se = se, mult = mult,
lower = lower, upper = upper, covers = covers, seed = seed(),
key = paste(n, level, input$mult))
})
# running total across "New sample" clicks; resets when n, level or multiplier change
running <- reactiveVal(list(key = "", covered = 0, total = 0))
observeEvent(sim(), {
s <- sim()
r <- running()
if (!identical(r$key, s$key)) r <- list(key = s$key, covered = 0, total = 0)
r$covered <- r$covered + sum(s$covers)
r$total <- r$total + 100
running(r)
})
output$ciPlot <- renderPlot({
s <- sim()
k <- sum(s$covers)
par(mar = c(4.5, 4.5, 3, 1))
plot(NA, xlim = range(c(s$lower, s$upper)), ylim = c(1, 100),
xlab = "Plant height (cm)", ylab = "Experiment",
main = paste(k, "of 100 intervals contain the true mean"))
segments(s$lower, 1:100, s$upper, 1:100,
col = ifelse(s$covers, "#002855", "#FF0000"),
lwd = ifelse(s$covers, 1.5, 3))
points(s$xbar, 1:100, pch = 16, cex = 0.6, col = "#002855")
abline(v = mu, col = "orange", lwd = 2)
})
output$stats <- renderPrint({
s <- sim()
r <- running()
lvl <- s$level
cat(sprintf("Intervals that contain 82: %d of 100\n", sum(s$covers)))
cat(sprintf("Mean interval width: %.2f cm\n\n", mean(s$upper - s$lower)))
cat(sprintf("Multiplier: qt(%.3f, %d) = %.3f vs qnorm(%.3f) = %.3f\n",
1 - s$a / 2, s$n - 1, qt(1 - s$a / 2, s$n - 1),
1 - s$a / 2, qnorm(1 - s$a / 2)))
cat(sprintf("Using: %s\n\n", if (input$mult == "t") "t" else "z"))
if (r$total > 0)
cat(sprintf("So far this session: %s of %s = %.1f%% (target %g%%)\n\n",
format(r$covered, big.mark = ","), format(r$total, big.mark = ","),
100 * r$covered / r$total, 100 * lvl))
ci1 <- t.test(s$m[1, ], conf.level = lvl)$conf.int
cat(sprintf("Experiment 1, by hand: %.2f to %.2f\n", s$lower[1], s$upper[1]))
cat(sprintf("Experiment 1, t.test(): %.2f to %.2f\n", ci1[1], ci1[2]))
if (input$mult == "t") cat("Same interval: t.test() uses the t multiplier.\n")
else cat("Different: t.test() uses qt(), not qnorm().\n")
})
output$code_out <- renderPrint({
s <- sim()
multline <- if (input$mult == "t")
sprintf("mult <- qt(1 - a / 2, df = n - 1) # %.3f", s$mult)
else
sprintf("mult <- qnorm(1 - a / 2) # %.3f", s$mult)
cat("# R code to reproduce this view\n")
cat("# (the running total across clicks is kept by the app only; it is not reproduced here)\n\n")
cat(sprintf("set.seed(%d)\n", s$seed))
cat(sprintf("n <- %d\n", s$n))
cat(sprintf("a <- 1 - %.2f # confidence level %g%%\n", s$level, 100 * s$level))
cat("m <- matrix(rnorm(100 * n, mean = 82, sd = 16), nrow = 100) # 100 experiments, one per row\n")
cat("xbar <- rowMeans(m)\n")
cat("se <- apply(m, 1, sd) / sqrt(n)\n")
cat(multline, "\n")
cat("lower <- xbar - mult * se\n")
cat("upper <- xbar + mult * se\n")
cat("covers <- lower <= 82 & upper >= 82\n")
cat("sum(covers) # intervals that contain the true mean\n")
cat("mean(upper - lower) # mean interval width\n\n")
cat("c(lower[1], upper[1]) # experiment 1 by hand\n")
cat(sprintf("t.test(m[1, ], conf.level = %.2f)$conf.int\n\n", s$level))
cat("plot(NA, xlim = range(c(lower, upper)), ylim = c(1, 100),\n")
cat(" xlab = \"Plant height (cm)\", ylab = \"Experiment\",\n")
cat(" main = paste(sum(covers), \"of 100 intervals contain the true mean\"))\n")
cat("segments(lower, 1:100, upper, 1:100,\n")
cat(" col = ifelse(covers, \"#002855\", \"#FF0000\"),\n")
cat(" lwd = ifelse(covers, 1.5, 3))\n")
cat("points(xbar, 1:100, pch = 16, cex = 0.6, col = \"#002855\")\n")
cat("abline(v = 82, col = \"orange\", lwd = 2)\n")
})
}
shinyApp(ui, server)t-test Explorer
Adjust the two group means, shared SD, and sample size. Both a Welch t-test and a Wilcoxon rank-sum test run on every sample. Toggle ranked view to see what the Wilcoxon test compares.
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 760
library(shiny)
ui <- fluidPage(
sidebarLayout(
sidebarPanel(
sliderInput("mean1", "Group A mean", min = 0, max = 20, value = 5, step = 0.5),
sliderInput("mean2", "Group B mean", min = 0, max = 20, value = 7, step = 0.5),
sliderInput("sd", "Std dev (both groups)", min = 0.5, max = 8, value = 2, step = 0.25),
sliderInput("n", "Sample size per group", min = 5, max = 200, value = 20, step = 5),
actionButton("resample", "New sample"),
hr(),
checkboxInput("show_ranks", "Show ranked values", value = FALSE),
checkboxInput("show_code", "Show R code", FALSE),
hr(),
helpText("Increase n with fixed means: the p-value drops even though the biology has not changed."),
helpText("Toggle ranked view: the Wilcoxon test compares rank sums, not raw means.")
),
mainPanel(
plotOutput("tPlot", height = "280px"),
verbatimTextOutput("tResult"),
conditionalPanel("input.show_code == true",
verbatimTextOutput("code_out")
)
)
)
)
server <- function(input, output, session) {
seed <- reactive(206 + input$resample)
dat <- reactive({
set.seed(seed())
a <- rnorm(input$n, input$mean1, input$sd)
b <- rnorm(input$n, input$mean2, input$sd)
list(a = a, b = b)
})
output$tPlot <- renderPlot({
d <- dat()
if (input$show_ranks) {
all_ranks <- rank(c(d$a, d$b))
plot_a <- all_ranks[seq_len(input$n)]
plot_b <- all_ranks[seq_len(input$n) + input$n]
ylab_text <- "Rank (pooled)"
main_text <- "Ranked values: what the Wilcoxon test compares"
} else {
plot_a <- d$a
plot_b <- d$b
ylab_text <- "Value"
main_text <- "Raw values: what the t-test compares"
}
ylo <- min(plot_a, plot_b) - 0.5
yhi <- max(plot_a, plot_b) + 0.5
par(mar = c(4, 4, 2, 1))
boxplot(list("Group A" = plot_a, "Group B" = plot_b),
col = c(adjustcolor("steelblue", 0.6),
adjustcolor("orange", 0.6)),
ylim = c(ylo, yhi),
ylab = ylab_text, main = main_text,
outline = FALSE, cex.axis = 1.2, cex.lab = 1.2)
stripchart(list("Group A" = plot_a, "Group B" = plot_b),
vertical = TRUE, method = "jitter",
jitter = 0.08, pch = 16, add = TRUE,
col = c(adjustcolor("steelblue", 0.5),
adjustcolor("orange", 0.5)))
points(1:2, c(mean(plot_a), mean(plot_b)),
pch = 18, cex = 2, col = "black")
})
output$tResult <- renderPrint({
d <- dat()
tres <- t.test(d$a, d$b)
wres <- wilcox.test(d$a, d$b, conf.int = TRUE, exact = FALSE)
sp <- sqrt((var(d$a) + var(d$b)) / 2) # pooled SD (equal n)
eff <- abs(mean(d$a) - mean(d$b)) / sp
cat(sprintf("Seed: %d\n\n", seed()))
cat("--- Welch t-test (parametric) ---\n")
cat(sprintf("t = %.3f, df = %.1f, p = %.4f %s\n",
tres$statistic, tres$parameter, tres$p.value,
ifelse(tres$p.value < 0.05, "[*]", "[ ]")))
cat(sprintf("95%% CI of difference: [%.3f, %.3f]\n",
tres$conf.int[1], tres$conf.int[2]))
cat(sprintf("Cohen's d (pooled SD): %.3f (%s)\n\n", eff,
ifelse(eff < 0.2, "negligible",
ifelse(eff < 0.5, "small",
ifelse(eff < 0.8, "medium", "large")))))
cat("--- Wilcoxon rank-sum test (non-parametric) ---\n")
cat(sprintf("W = %.0f, p = %.4f %s\n",
wres$statistic, wres$p.value,
ifelse(wres$p.value < 0.05, "[*]", "[ ]")))
cat(sprintf("95%% CI of location shift: [%.3f, %.3f]\n",
wres$conf.int[1], wres$conf.int[2]))
cat(sprintf("\nAgree on significance: %s\n",
ifelse((tres$p.value < 0.05) == (wres$p.value < 0.05),
"YES", "NO, they disagree")))
})
output$code_out <- renderPrint({
cat("# R code to reproduce this view:\n\n")
cat(sprintf("set.seed(%d)\n", seed()))
cat(sprintf("a <- rnorm(%d, mean = %g, sd = %g) # group A\n",
input$n, input$mean1, input$sd))
cat(sprintf("b <- rnorm(%d, mean = %g, sd = %g) # group B\n\n",
input$n, input$mean2, input$sd))
cat("t.test(a, b) # Welch t-test\n")
cat("wilcox.test(a, b, conf.int = TRUE, exact = FALSE) # Wilcoxon rank-sum\n")
cat("sp <- sqrt((var(a) + var(b)) / 2) # pooled SD (equal n)\n")
cat("abs(mean(a) - mean(b)) / sp # Cohen's d\n\n")
if (input$show_ranks) {
cat("# ranked view: pool both groups and rank them\n")
cat("r <- rank(c(a, b))\n")
cat(sprintf("a <- r[1:%d]; b <- r[%d:%d]\n\n", input$n, input$n + 1, 2 * input$n))
}
cat("boxplot(list(\"Group A\" = a, \"Group B\" = b), outline = FALSE,\n")
cat(" col = c(\"steelblue\", \"orange\"))\n")
cat("stripchart(list(\"Group A\" = a, \"Group B\" = b), vertical = TRUE,\n")
cat(" method = \"jitter\", pch = 16, add = TRUE)\n")
cat("points(1:2, c(mean(a), mean(b)), pch = 18, cex = 2)\n")
})
}
shinyApp(ui, server)p-values, Power, and Multiple Testing
This is the watered vs droughted wheat trial from lecture, run 500 times. Start with no drought effect and look at the spread of p-values: about 5% land below 0.05 by chance. Then add an effect and watch power grow with n, and in the second tab see why many tests need a correction (Week 1 script, sections 13, 16 and 26).
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 900
library(shiny)
SD <- 12 # plant-to-plant SD (cm), fixed as in lecture
MU <- 84 # watered mean height (cm)
NSIM <- 500 # experiments in tab 1
ui <- fluidPage(
sidebarLayout(
sidebarPanel(
sliderInput("effect", "True drought effect (cm)", 0, 20, 0, 1),
sliderInput("n", "Plants per group (n)", 3, 60, 15, 1),
helpText("Watered plants ~ N(84, 12), droughted plants ~ N(84 - effect, 12). The SD is fixed at 12 cm."),
actionButton("resample", "New sample"),
hr(),
checkboxInput("show_code", "Show R code", FALSE),
helpText("Try: effect = 0, then 10 cm with n = 15 (the lecture's power example), then n = 30.")
),
mainPanel(
tabsetPanel(id = "tab",
tabPanel("One test, 500 experiments", value = "one",
plotOutput("histPlot", height = "380px"),
verbatimTextOutput("oneStats")
),
tabPanel("Many tests, one study", value = "many",
fluidRow(
column(6, sliderInput("m", "Number of tests", 5, 100, 20, 1)),
column(6, sliderInput("k", "How many are real effects", 0, 10, 2, 1))
),
helpText("Real tests use the drought effect slider. If it is set to 0, they use 12 cm."),
plotOutput("manyPlot", height = "430px"),
verbatimTextOutput("manyStats")
)
),
conditionalPanel("input.show_code == true",
verbatimTextOutput("code_out")
)
)
)
)
server <- function(input, output, session) {
seed <- reactive(206 + input$resample)
# ---- tab 1: 500 experiments, Welch t by hand (vectorised) ----
sims <- reactive({
n <- input$n; effect <- input$effect
set.seed(seed())
w <- matrix(rnorm(NSIM * n, MU, SD), NSIM)
d <- matrix(rnorm(NSIM * n, MU - effect, SD), NSIM)
vw <- apply(w, 1, var); vd <- apply(d, 1, var)
se <- sqrt(vw / n + vd / n)
t <- (rowMeans(w) - rowMeans(d)) / se
df <- se^4 / ((vw / n)^2 / (n - 1) + (vd / n)^2 / (n - 1))
p <- 2 * pt(-abs(t), df)
list(p = p, check = t.test(w[1, ], d[1, ])$p.value)
})
output$histPlot <- renderPlot({
p <- sims()$p
hit_col <- if (input$effect == 0) "#FF0000" else "#00B050"
cols <- c(hit_col, rep(adjustcolor("steelblue", 0.7), 19))
par(mar = c(4.5, 4.5, 3, 1))
h <- hist(p, breaks = seq(0, 1, 0.05), plot = FALSE)
plot(h, col = cols, border = "white",
ylim = c(0, max(h$counts, 25) * 1.1),
xlab = "p-value", ylab = "Number of experiments",
main = sprintf("p-values from %d experiments (effect = %g cm, n = %d)",
NSIM, input$effect, input$n))
abline(h = NSIM * 0.05, lty = 2, lwd = 2)
legend("topright", bty = "n",
legend = c(if (input$effect == 0) "p < 0.05: false positives"
else "p < 0.05: true positives",
"Expected count per bar if no effect (25)"),
fill = c(hit_col, NA), border = c("white", NA),
lty = c(NA, 2), lwd = c(NA, 2))
})
output$oneStats <- renderPrint({
s <- sims()
cat(sprintf("Share of experiments with p < 0.05: %.3f (%d of %d)\n",
mean(s$p < 0.05), sum(s$p < 0.05), NSIM))
if (input$effect == 0) {
cat("Effect = 0: expect 0.05, the Type I error rate.\n")
} else {
pw <- power.t.test(n = input$n, delta = input$effect, sd = SD)$power
cat(sprintf("power.t.test(n = %d, delta = %g, sd = 12)$power = %.3f\n",
input$n, input$effect, pw))
cat(sprintf("Type II error rate (missed effects): %.3f\n", mean(s$p >= 0.05)))
}
})
# ---- tab 2: m tests in one study, k of them real ----
many <- reactive({
n <- input$n; m <- input$m; k <- min(input$k, m)
eff <- if (input$effect == 0) 12 else input$effect
set.seed(seed())
p_null <- if (m - k > 0) replicate(m - k, t.test(rnorm(n, MU, SD), rnorm(n, MU, SD))$p.value)
p_real <- if (k > 0) replicate(k, t.test(rnorm(n, MU, SD), rnorm(n, MU - eff, SD))$p.value)
p <- c(p_null, p_real)
real <- c(rep(FALSE, m - k), rep(TRUE, k))
list(p = p, real = real, m = m, k = k, eff = eff,
raw = p < 0.05,
holm = p.adjust(p, "holm") < 0.05,
bh = p.adjust(p, "BH") < 0.05)
})
output$manyPlot <- renderPlot({
r <- many(); m <- r$m
pp <- pmax(r$p, 1e-300)
col <- ifelse(r$holm, "#7030A0",
ifelse(r$bh, "#0070C0",
ifelse(r$raw, "#FF0000", "grey60")))
y <- seq_len(m)
par(mar = c(4.5, 4.5, 3, 1))
plot(pp, y, log = "x", type = "n",
xlim = c(min(pp, 1e-4), 1), ylim = c(0.5, m + 0.5),
yaxt = "n", xlab = "p-value (log scale)", ylab = "Test",
main = sprintf("%d tests, %d real (effect %g cm), n = %d per group",
m, r$k, r$eff, input$n))
axis(2, at = pretty(y), las = 1)
abline(h = y, col = "grey92")
abline(v = 0.05, lty = 2, lwd = 2, col = "#FF0000")
abline(v = 0.05 / m, lty = 2, lwd = 2, col = "#7030A0")
points(pp, y, pch = 16, cex = 1.6, col = col)
# mark which tests are real effects
points(rep(1, m), y, pch = ifelse(r$real, 17, 1),
col = ifelse(r$real, "#00B050", "grey70"), cex = 1.2)
legend("bottomleft", bty = "n", cex = 0.85,
legend = c("significant after Holm", "after BH (not Holm)",
"raw p < 0.05 only", "not significant",
"real effect (right edge)", "no effect",
"p = 0.05", sprintf("Bonferroni cut 0.05/%d", m)),
pch = c(16, 16, 16, 16, 17, 1, NA, NA),
lty = c(NA, NA, NA, NA, NA, NA, 2, 2),
col = c("#7030A0", "#0070C0", "#FF0000", "grey60",
"#00B050", "grey70", "#FF0000", "#7030A0"))
})
output$manyStats <- renderPrint({
r <- many()
tab <- rbind(
"raw p < 0.05" = c(sum(r$raw & r$real), sum(r$raw & !r$real)),
"Holm" = c(sum(r$holm & r$real), sum(r$holm & !r$real)),
"BH (FDR)" = c(sum(r$bh & r$real), sum(r$bh & !r$real)))
colnames(tab) <- c("true hits", "false hits")
print(tab)
if (r$k == 0) cat("\nNo real effects: every hit is a false positive.\n")
else cat(sprintf("\nReal effects: %d of %d tests (the last %d).\n", r$k, r$m, r$k))
cat(sprintf("Chance of at least one false positive if all %d tests were null:\n", r$m))
cat(sprintf("1 - 0.95^%d = %.3f\n", r$m, 1 - 0.95^r$m))
})
output$code_out <- renderPrint({
n <- input$n; s <- seed()
if (input$tab == "one") {
effect <- input$effect
cat("# Reproduce this view: 500 experiments, watered vs droughted\n\n")
cat(sprintf("set.seed(%d)\n", s))
cat(sprintf("n <- %d; effect <- %g\n", n, effect))
cat("w <- matrix(rnorm(500 * n, 84, 12), 500) # watered, one row per experiment\n")
cat("d <- matrix(rnorm(500 * n, 84 - effect, 12), 500) # droughted\n\n")
cat("# Welch t and p by hand for all 500 rows at once\n")
cat("# (same as t.test(w[i, ], d[i, ])$p.value, just faster)\n")
cat("vw <- apply(w, 1, var); vd <- apply(d, 1, var)\n")
cat("se <- sqrt(vw / n + vd / n)\n")
cat("t <- (rowMeans(w) - rowMeans(d)) / se\n")
cat("df <- se^4 / ((vw / n)^2 / (n - 1) + (vd / n)^2 / (n - 1))\n")
cat("p <- 2 * pt(-abs(t), df)\n\n")
cat("mean(p < 0.05) # share significant\n")
cat("hist(p, breaks = seq(0, 1, 0.05))\n")
cat("abline(h = 500 * 0.05, lty = 2) # expected per bar if no effect\n\n")
if (effect > 0) {
cat("power.t.test(n = n, delta = effect, sd = 12)$power\n\n")
} else {
cat("# effect = 0: expect mean(p < 0.05) near 0.05, the Type I error rate\n\n")
}
cat("# check: the by-hand p equals t.test()\n")
cat("all.equal(t.test(w[1, ], d[1, ])$p.value, p[1])\n")
} else {
r <- many()
cat("# Reproduce this view: many tests in one study\n\n")
cat(sprintf("set.seed(%d)\n", s))
cat(sprintf("n <- %d; effect <- %g%s\n", n, r$eff,
if (input$effect == 0) " # slider at 0, so real tests use 12 cm" else ""))
null_line <- sprintf("replicate(%d, t.test(rnorm(n, 84, 12), rnorm(n, 84, 12))$p.value)",
r$m - r$k)
real_line <- sprintf("replicate(%d, t.test(rnorm(n, 84, 12), rnorm(n, 84 - effect, 12))$p.value)",
r$k)
if (r$k == 0) {
cat("p <-", null_line, "\n")
} else if (r$k == r$m) {
cat("p <-", real_line, "\n")
} else {
cat("p <- c(", null_line, ",\n ", real_line, ")\n", sep = "")
}
if (r$k == r$m)
cat(sprintf("# all %d tests are real effects\n\n", r$m))
else if (r$k > 0)
cat(sprintf("# the last %d tests (%d to %d) are real effects\n\n",
r$k, r$m - r$k + 1, r$m))
else cat("# no real effects: every hit is false\n\n")
cat("which(p < 0.05) # raw\n")
cat("which(p.adjust(p, \"holm\") < 0.05) # after Holm\n")
cat("which(p.adjust(p, \"BH\") < 0.05) # after BH (FDR)\n")
cat(sprintf("1 - 0.95^%d # P(at least one false positive)\n\n",
r$m))
cat("plot(p, seq_along(p), log = \"x\", pch = 16, xlab = \"p-value\", ylab = \"Test\")\n")
cat(sprintf("abline(v = c(0.05, 0.05 / %d), lty = 2) # raw and Bonferroni cuts\n", r$m))
}
})
}
shinyApp(ui, server)Correlation Explorer
Set the true correlation strength, sample size, and number of outliers. Watch how Pearson gets pulled by outliers more than Spearman does: Spearman works on ranks, so it sees only the order of the points, not how far out an outlier sits.
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 640
library(shiny)
ui <- fluidPage(
sidebarLayout(
sidebarPanel(
sliderInput("r", "True correlation (r)", min = -0.99, max = 0.99, value = 0.5, step = 0.05),
sliderInput("n", "Sample size (n)", min = 10, max = 300, value = 50, step = 10),
sliderInput("out", "Outliers", min = 0, max = 5, value = 0, step = 1),
actionButton("resample", "New sample"),
hr(),
checkboxInput("show_code", "Show R code", FALSE),
helpText("Add outliers to see why Spearman is more robust than Pearson.")
),
mainPanel(
plotOutput("corrPlot", height = "300px"),
verbatimTextOutput("corrResult"),
conditionalPanel("input.show_code == true",
verbatimTextOutput("code_out")
)
)
)
)
server <- function(input, output, session) {
seed <- reactive(206 + input$resample)
dat <- reactive({
set.seed(seed())
r_val <- max(-0.99, min(0.99, input$r))
sigma <- matrix(c(1, r_val, r_val, 1), 2, 2)
L <- chol(sigma)
raw <- matrix(rnorm(input$n * 2), ncol = 2) %*% L
df <- data.frame(x = raw[,1], y = raw[,2], outlier = FALSE)
if (input$out > 0) {
outs <- data.frame(x = runif(input$out, 2, 4),
y = runif(input$out, -4, -2),
outlier = TRUE)
df <- rbind(df, outs)
}
df
})
output$corrPlot <- renderPlot({
d <- dat()
norm <- d[!d$outlier, ]
par(mar = c(4, 4, 1, 1))
plot(d$x, d$y,
col = ifelse(d$outlier, "red", adjustcolor("gray40", 0.7)),
pch = 16, cex = 1.3,
xlab = "X", ylab = "Y",
cex.axis = 1.2, cex.lab = 1.2)
if (nrow(norm) > 1)
abline(lm(y ~ x, data = norm), col = "steelblue", lwd = 2)
abline(lm(y ~ x, data = d), col = "red", lty = 2, lwd = 1.8)
legend("topleft",
legend = c("Normal pts", "Outlier", "Fit (no outliers)", "Fit (all)"),
col = c("gray40", "red", "steelblue", "red"),
pch = c(16, 16, NA, NA),
lty = c(NA, NA, 1, 2),
lwd = c(NA, NA, 2, 1.8),
bty = "n", cex = 0.9)
})
output$corrResult <- renderPrint({
d <- dat()
pear <- cor.test(d$x, d$y, method = "pearson")
spear <- cor.test(d$x, d$y, method = "spearman", exact = FALSE)
cat(sprintf("Pearson r = %+.3f (p = %.4f)\n", pear$estimate, pear$p.value))
cat(sprintf("Spearman ρ = %+.3f (p = %.4f)\n", spear$estimate, spear$p.value))
cat(sprintf("\nTrue r = %.2f | n = %d | outliers = %d | seed = %d\n",
input$r, input$n, input$out, seed()))
})
output$code_out <- renderPrint({
r_val <- max(-0.99, min(0.99, input$r))
cat("# R code to reproduce this view:\n\n")
cat(sprintf("set.seed(%d)\n", seed()))
cat(sprintf("r <- %g # true correlation\n", r_val))
cat("L <- chol(matrix(c(1, r, r, 1), 2, 2))\n")
cat(sprintf("raw <- matrix(rnorm(%d * 2), ncol = 2) %%*%% L # correlated pairs\n",
input$n))
cat("x <- raw[, 1]\ny <- raw[, 2]\n")
cat("x_clean <- x; y_clean <- y # before outliers\n")
if (input$out > 0) {
cat(sprintf("x <- c(x, runif(%d, 2, 4)) # outliers: high x ...\n", input$out))
cat(sprintf("y <- c(y, runif(%d, -4, -2)) # ... and low y\n", input$out))
}
cat("\ncor.test(x, y, method = \"pearson\")\n")
cat("cor.test(x, y, method = \"spearman\", exact = FALSE)\n\n")
cat("plot(x, y, pch = 16)\n")
cat("abline(lm(y_clean ~ x_clean), col = \"steelblue\", lwd = 2) # fit without outliers\n")
cat("abline(lm(y ~ x), col = \"red\", lty = 2, lwd = 2) # fit with all points\n")
})
}
shinyApp(ui, server)